# -*- coding: utf-8 -*-
# coding=utf-8
import os
from osgeo import gdal, ogr, osr
import math
import numpy as np
import json
import codecs

'''把栅格转为矢量面
    # extend为获取范围，默认参数是北京的矩形范围
    isBj为标识是否只计算北京范围
    inrasterPath为输入的tiff路径
    shpOutPath为输出的shp路径
'''
inrasterPath='/data/py_tiff_back/testData/tiffData/RS_FOR_BJ_FMP_202409030610_202409031110_1H_5_1_1.tif'
shpOutPath='/data/py_tiff_back/testData/outData/test01.shp'
def createGeo(isBj=1):
    if (isBj != 1):
        extent = (0, 180, 0, 80)
    else:
        extent = (115.41730212, 117.50017715999999, 39.438288, 41.05923498)
    # inraster = gdal.Open('G:\\code\\ydh\\py_tiff_back\\testData\\tiffData\\clipped_image_0.tif')
    inraster = gdal.Open(inrasterPath)
    # prj = osr.SpatialReference()
    # prj.ImportFromWkt(inraster.GetProjection())
    #
    driver=ogr.GetDriverByName('ESRI Shapefile')
    ds=driver.CreateDataSource(shpOutPath)


    prj = osr.SpatialReference()
    prj.ImportFromWkt(
        'GEOGCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID["CGCS2000",6378137,298.257222101004,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4490"]]'
    )

    # 图层
    layer=ds.CreateLayer('test', srs=prj,geom_type=ogr.wkbPolygon)
    # 字段，图层把字段加上
    fieldDefn=ogr.FieldDefn('id',ogr.OFTString)
    fieldDefn.SetWidth(8)
    layer.CreateField(fieldDefn)
    fieldValue = ogr.FieldDefn('value', ogr.OFTReal)
    layer.CreateField(fieldValue)
    # 图层获取属性集合对象feature，方便layer后面往里面加feature
    featureDefn=layer.GetLayerDefn()
    feature=ogr.Feature(featureDefn)

    inband = inraster.GetRasterBand(1)
    inbandData = inband.ReadAsArray()
    A_geotrans = inraster.GetGeoTransform()

    for i in  range(0,inbandData.shape[0]):
        for j in range(0, inbandData.shape[1]):
            if((A_geotrans[0] + A_geotrans[1] * j) > extent[1] or (A_geotrans[0] + A_geotrans[1] * j) < extent[0]):
                continue
            elif((A_geotrans[3] + A_geotrans[5] * i) > extent[3] or (A_geotrans[3] + A_geotrans[5] * i) < extent[2]):
                continue
            elif(not np.isnan(inbandData[i][j]) and math.ceil(inbandData[i][j])>0):
                # print(str(i) + '_' + str(j))
                polygon = None
                outring = None
                polygon = ogr.Geometry(ogr.wkbPolygon)
                outring = ogr.Geometry(ogr.wkbLinearRing)

                Lon = A_geotrans[0] + A_geotrans[1] * j
                Lat = A_geotrans[3] + A_geotrans[5] * i

                # outring.AddPoint(120, 30)
                # outring.AddPoint(121, 30)
                # outring.AddPoint(121, 29)
                # outring.AddPoint(120, 29)
                # outring.AddPoint(120, 30)
                outring.AddPoint(Lon, Lat)
                outring.AddPoint(Lon + A_geotrans[1], Lat)
                outring.AddPoint(Lon + A_geotrans[1], Lat + A_geotrans[5])
                outring.AddPoint(Lon, Lat + A_geotrans[5])
                outring.AddPoint(Lon, Lat)
                polygon.AddGeometry(outring)
                feature.SetGeometry(polygon)
                feature.SetField('id', str(i) + '_' + str(j))
                feature.SetField('value', round(float(inbandData[i][j]),2))
                layer.CreateFeature(feature)


    ds.Destroy()



createGeo()


def openShp():
    inshp="G:\\code\\ydh\\py_tiff_back\\testData\\shp\\beijing4490.shp"
    dataSource=ogr.Open(inshp)
    layer=dataSource.GetLayer(0)
    spPro=layer.GetSpatialRef()
    extent=layer.GetExtent()
    print(spPro)
    print(extent)
    print(extent[0])
    print(extent[1])
# openShp()

def openTiff(path):
    inraster = gdal.Open(path)  # 读取路径中的栅格数据
    # prj = osr.SpatialReference()
    # prj.ImportFromWkt(inraster.GetProjection())
    prj = inraster.GetProjection() # str类型的proj
    A_geotrans = inraster.GetGeoTransform()


    # inband.CreateMaskBand(0)
    prj = osr.SpatialReference()
    prj.ImportFromWkt(inraster.GetProjection())  # 读取栅格数据的投影信息，用来为后面生成的矢量做准备
    #
    outshp = path[:-4] + "5.shp"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outshp):  # 若文件已经存在，则删除它继续重新做一遍
        drv.DeleteDataSource(outshp)
    Polygon = drv.CreateDataSource(outshp)  # 创建一个目标文件
    Poly_layer = Polygon.CreateLayer(path, srs=prj, geom_type=ogr.wkbMultiPolygon)
    newField = ogr.FieldDefn('value',
                             ogr.OFTInteger)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value,OFTReal浮点型，ogr.OFTInteger是整数型
    Poly_layer.CreateField(newField)

def tiffDo(tiffName):

    return None;
